I need your help!

I want your feedback to make the book better for you and other readers. If you find typos, errors, or places where the text may be improved, please let me know. The best ways to provide feedback are by GitHub or hypothes.is annotations.

Opening an issue or submitting a pull request on GitHub: https://github.com/isaactpetersen/Fantasy-Football-Analytics-Textbook

Adding an annotation using hypothes.is. To add an annotation, select some text and then click the symbol on the pop-up menu. To see the annotations of others, click the symbol in the upper right-hand corner of the page.

11  Mixed Models

11.1 Getting Started

11.1.1 Load Packages

library("lme4")
library("lmerTest")
library("MuMIn")
library("emmeans")
library("sjstats")
library("mgcv")
library("AICcmodavg")
library("bbmle")
library("parallel")
library("plotly")
library("viridis")
library("tidyverse")

11.1.2 Specify Package Options

emm_options(lmerTest.limit = 100000)
emm_options(pbkrtest.limit = 100000)

11.1.3 Load Data

load(file = "./data/nfl_depthCharts.RData")
load(file = "./data/player_stats_weekly.RData")
load(file = "./data/player_stats_seasonal.RData")

We created the player_stats_weekly.RData and player_stats_seasonal.RData objects in Section 3.21.3.

11.2 Overview of Mixed Models

We will discuss a modeling framework that goes by many terms, including mixed models, mixed-effects models, multilevel models, hierarchical linear models.

11.2.1 Ecological Fallacy

11.2.2 Simpson’s Paradox

11.3 Fantasy Points Per Season by Position, Age, and Experience

player_stats_seasonal_offense_subset <- player_stats_seasonal_offense %>% 
  filter(position_group %in% c("QB","RB","WR","TE"))

player_stats_seasonal_offense_subset$position[which(player_stats_seasonal_offense_subset$position == "HB")] <- "RB"

player_stats_seasonal_kicking_subset <- player_stats_seasonal_kicking %>% 
  filter(position == "K")

player_stats_seasonal_offense_subset <- bind_rows(
  player_stats_seasonal_offense_subset,
  player_stats_seasonal_kicking_subset
)

player_stats_seasonal_offense_subset$player_idFactor <- factor(player_stats_seasonal_offense_subset$player_id)
player_stats_seasonal_offense_subset$positionFactor <- factor(player_stats_seasonal_offense_subset$position)

seasons17week <- 2001:2020
seasons18week <- 2021:max(nfl_depthCharts$season, na.rm = TRUE)

endOfSeasonDepthCharts <- nfl_depthCharts %>% 
  filter((season %in% seasons17week & week == 18) | (season %in% seasons18week & week == 19)) # get end-of-season depth charts

qb1s <- endOfSeasonDepthCharts %>% 
  filter(position == "QB", depth_team == 1)

fb1s <- endOfSeasonDepthCharts %>% 
  filter(position == "FB", depth_team == 1)

k1s <- endOfSeasonDepthCharts %>% 
  filter(position == "K", depth_team == 1)

rb1s <- endOfSeasonDepthCharts %>% 
  filter(position == "RB", depth_team == 1)

wr1s <- endOfSeasonDepthCharts %>% 
  filter(position == "WR", depth_team == 1)

te1s <- endOfSeasonDepthCharts %>% 
  filter(position == "TE", depth_team == 1)

player_stats_seasonal_offense_subsetDepth <- player_stats_seasonal_offense_subset %>% 
  filter(player_id %in% c(
    qb1s$gsis_id,
    fb1s$gsis_id,
    k1s$gsis_id,
    rb1s$gsis_id,
    wr1s$gsis_id,
    te1s$gsis_id
    ))

Create a newdata object for generating the plots of model-implied fantasy points by age and position:

pointsPerSeason_positionAge_newData <- expand.grid(
  positionFactor = factor(c("FB","QB","RB","TE","WR")), #,"K"
  age = seq(from = 20, to = 40, length.out = 10000)
)

pointsPerSeason_positionAge_newData$ageCentered20 <- pointsPerSeason_positionAge_newData$age - 20
pointsPerSeason_positionAge_newData$ageCentered20Quadratic <- pointsPerSeason_positionAge_newData$ageCentered20 ^ 2
pointsPerSeason_positionAge_newData$years_of_experience <- floor(pointsPerSeason_positionAge_newData$age - 22) # assuming that most players start at age 22 (i.e., rookie year) and thus have 1 year of experience at age 23
pointsPerSeason_positionAge_newData$years_of_experience[which(pointsPerSeason_positionAge_newData$years_of_experience < 0)] <- 0

Create an object with complete cases for generating the plots of individuals’ -implied fantasy points by age and position:

player_stats_seasonal_offense_subsetCC <- player_stats_seasonal_offense_subset %>%
  filter(
    !is.na(player_idFactor),
    !is.na(fantasy_points),
    !is.na(positionFactor),
    !is.na(ageCentered20),
    !is.na(ageCentered20Quadratic),
    !is.na(years_of_experience))

11.3.1 Scatterplots of Fantasy Points by Age and Position

Scatterplots are a helpful tool for quickly examining the association between two variables. However, scatterplots—as well as correlation and multiple regression—can hide meaningful associations that differ across units of analysis.

11.3.1.1 Quarterbacks

A scatterplot of Quarterbacks’ fantasy points by age is in Figure 11.1.

plot_scatterplotFantasyPointsByAgeQB <- ggplot(
  data = player_stats_seasonal_offense_subset %>% 
    filter(position == "QB") %>% 
    mutate(
      age = round(age, 2),
      fantasy_points = round(fantasy_points, 2)
    ),
  mapping = aes(
    x = age,
    y = fantasy_points,
    color = player_id)) +
  geom_point(
    aes(
      text = player_display_name, # add player name for mouse over tooltip
      label = season # add season for mouse over tooltip
  )) +
  geom_smooth(
    mapping = aes(
    x = age,
    y = fantasy_points),
    inherit.aes = FALSE
  ) +
  scale_color_viridis(discrete = TRUE) +
  labs(
    x = "Player Age (years)",
    y = "Fantasy Points (Season)",
    title = "Fantasy Points (Season) by Player Age: Quarterbacks"
  ) +
  theme_classic() +
  theme(legend.position = "none")

ggplotly(
  plot_scatterplotFantasyPointsByAgeQB,
  tooltip = c("age","fantasy_points","text","label"))

Figure 11.1: Scatterplot of Fantasy Points by Age for Quarterbacks.

Based on the scatterplot (and the bivariate association below), Quarterbacks’ fantasy points appear to increase with age.

cor.test(
  formula = ~ age + fantasy_points,
  data = player_stats_seasonal_offense_subset %>% filter(position == "QB")
)

    Pearson's product-moment correlation

data:  age and fantasy_points
t = 3.2358, df = 1051, p-value = 0.001251
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 0.03914186 0.15877861
sample estimates:
       cor 
0.09931915 

11.3.1.2 Fullbacks

A scatterplot of Fullbacks’ fantasy points by age is in Figure 11.2.

plot_scatterplotFantasyPointsByAgeFB <- ggplot(
  data = player_stats_seasonal_offense_subset %>% 
    filter(position == "FB") %>% 
    mutate(
      age = round(age, 2),
      fantasy_points = round(fantasy_points, 2)
    ),
  mapping = aes(
    x = age,
    y = fantasy_points,
    color = player_id)) +
  geom_point(
    aes(
      text = player_display_name, # add player name for mouse over tooltip
      label = season # add season for mouse over tooltip
  )) +
  geom_smooth(
    mapping = aes(
    x = age,
    y = fantasy_points),
    inherit.aes = FALSE
  ) +
  scale_color_viridis(discrete = TRUE) +
  labs(
    x = "Player Age (years)",
    y = "Fantasy Points (Season)",
    title = "Fantasy Points (Season) by Player Age: Fullbacks"
  ) +
  theme_classic() +
  theme(legend.position = "none")

ggplotly(
  plot_scatterplotFantasyPointsByAgeFB,
  tooltip = c("age","fantasy_points","text","label"))

Figure 11.2: Scatterplot of Fantasy Points by Age for Fullbacks.

Based on the scatterplot, Fullbacks’ fantasy points appear to be relatively stable across ages.

11.3.1.3 Running Backs

A scatterplot of Running Backs’ fantasy points by age is in Figure 11.3.

plot_scatterplotFantasyPointsByAgeRB <- ggplot(
  data = player_stats_seasonal_offense_subset %>% 
    filter(position == "RB") %>% 
    mutate(
      age = round(age, 2),
      fantasy_points = round(fantasy_points, 2)
    ),
  mapping = aes(
    x = age,
    y = fantasy_points,
    color = player_id)) +
  geom_point(
    aes(
      text = player_display_name, # add player name for mouse over tooltip
      label = season # add season for mouse over tooltip
  )) +
  geom_smooth(
    mapping = aes(
    x = age,
    y = fantasy_points),
    inherit.aes = FALSE
  ) +
  scale_color_viridis(discrete = TRUE) +
  labs(
    x = "Player Age (years)",
    y = "Fantasy Points (Season)",
    title = "Fantasy Points (Season) by Player Age: Running Backs"
  ) +
  theme_classic() +
  theme(legend.position = "none")

ggplotly(
  plot_scatterplotFantasyPointsByAgeRB,
  tooltip = c("age","fantasy_points","text","label"))

Figure 11.3: Scatterplot of Fantasy Points by Age for Running Backs.

Based on the scatterplot, Running Backs’ fantasy points appear to be relatively stable across ages.

11.3.1.4 Wide Receivers

A scatterplot of Wide Receivers’ fantasy points by age is in Figure 11.4.

plot_scatterplotFantasyPointsByAgeWR <- ggplot(
  data = player_stats_seasonal_offense_subset %>% 
    filter(position == "WR") %>% 
    mutate(
      age = round(age, 2),
      fantasy_points = round(fantasy_points, 2)
    ),
  mapping = aes(
    x = age,
    y = fantasy_points,
    color = player_id)) +
  geom_point(
    aes(
      text = player_display_name, # add player name for mouse over tooltip
      label = season # add season for mouse over tooltip
  )) +
  geom_smooth(
    mapping = aes(
    x = age,
    y = fantasy_points),
    inherit.aes = FALSE
  ) +
  scale_color_viridis(discrete = TRUE) +
  labs(
    x = "Player Age (years)",
    y = "Fantasy Points (Season)",
    title = "Fantasy Points (Season) by Player Age: Wide Receivers"
  ) +
  theme_classic() +
  theme(legend.position = "none")

ggplotly(
  plot_scatterplotFantasyPointsByAgeWR,
  tooltip = c("age","fantasy_points","text","label"))

Figure 11.4: Scatterplot of Fantasy Points by Age for Wide Receivers.

Based on the scatterplot, Wide Receivers’ fantasy points appear to be relatively stable across ages.

11.3.1.5 Tight Ends

A scatterplot of Tight Ends’ fantasy points by age is in Figure 11.5.

plot_scatterplotFantasyPointsByAgeTE <- ggplot(
  data = player_stats_seasonal_offense_subset %>% 
    filter(position == "TE") %>% 
    mutate(
      age = round(age, 2),
      fantasy_points = round(fantasy_points, 2)
    ),
  mapping = aes(
    x = age,
    y = fantasy_points,
    color = player_id)) +
  geom_point(
    aes(
      text = player_display_name, # add player name for mouse over tooltip
      label = season # add season for mouse over tooltip
  )) +
  geom_smooth(
    mapping = aes(
    x = age,
    y = fantasy_points),
    inherit.aes = FALSE
  ) +
  scale_color_viridis(discrete = TRUE) +
  labs(
    x = "Player Age (years)",
    y = "Fantasy Points (Season)",
    title = "Fantasy Points (Season) by Player Age: Tight Ends"
  ) +
  theme_classic() +
  theme(legend.position = "none")

ggplotly(
  plot_scatterplotFantasyPointsByAgeTE,
  tooltip = c("age","fantasy_points","text","label"))

Figure 11.5: Scatterplot of Fantasy Points by Age for Tight Ends.

Based on the scatterplot (and the bivariate association below), Tight Ends’ fantasy points appear to increase with age.

cor.test(
  formula = ~ age + fantasy_points,
  data = player_stats_seasonal_offense_subset %>% filter(position == "TE")
)

    Pearson's product-moment correlation

data:  age and fantasy_points
t = 5.3884, df = 1341, p-value = 8.388e-08
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 0.09280906 0.19753022
sample estimates:
      cor 
0.1455774 

11.3.2 Plots of Raw Trajectories of Fantasy Points By Age and Player

Scatterplots can be helpful for quickly visualizing the association between two variables. However, as mentioned earlier, scatterplots can hide the association between variables at different units of analysis. For instance, consider that we are trying to predict how a player will perform based on their age. We are interested not only in what the association is between age and fantasy points between players (i.e., a between-person association). We are also interested in what the association is between age and fantasy points within a given player (and within each player; i.e., a within-individual association). Arguably, the within-individual association between age and fantasy points is more relevant to the prediction of performance than the association between age and fantasy points between players. Assuming that the between-player association between age and fantasy points is the same as the within-player association when it is not is an example of the ecological fallacy.

Below, we depict players’ raw trajectories of fantasy points as a function of age. These are known as spaghetti plots. By examining the trajectory for each player, we can get a better understanding of hor performance changes (within an individual) as a function of age.

11.3.2.1 Quarterbacks

A plot of Quarterbacks’ raw fantasy points data by age is in Figure 11.6.

plot_rawFantasyPointsByAgeQB <- ggplot(
  data = player_stats_seasonal_offense_subset %>% 
    filter(position == "QB") %>% 
    mutate(
      age = round(age, 2),
      fantasy_points = round(fantasy_points, 2)
    ),
  mapping = aes(
    x = age,
    y = fantasy_points,
    color = player_id)) +
  geom_line(
    aes(
      text = player_display_name, # add player name for mouse over tooltip
      label = season # add season for mouse over tooltip
  )) +
  scale_color_viridis(discrete = TRUE) +
  labs(
    x = "Player Age (years)",
    y = "Fantasy Points (Season)",
    title = "Fantasy Points (Season) by Player Age: Quarterbacks"
  ) +
  theme_classic() +
  theme(legend.position = "none")

ggplotly(
  plot_rawFantasyPointsByAgeQB,
  tooltip = c("age","fantasy_points","text","label"))

Figure 11.6: Plot of Raw Trajectories of Fantasy Points by Age for Quarterbacks.

11.3.2.2 Fullbacks

A plot of Fullbacks’ raw fantasy points data by age is in Figure 11.7.

plot_rawFantasyPointsByAgeFB <- ggplot(
  data = player_stats_seasonal_offense_subset %>% 
    filter(position == "FB") %>% 
    mutate(
      age = round(age, 2),
      fantasy_points = round(fantasy_points, 2)
    ),
  mapping = aes(
    x = age,
    y = fantasy_points,
    color = player_id)) +
  geom_line(
    aes(
      text = player_display_name, # add player name for mouse over tooltip
      label = season # add season for mouse over tooltip
  )) +
  scale_color_viridis(discrete = TRUE) +
  labs(
    x = "Player Age (years)",
    y = "Fantasy Points (Season)",
    title = "Fantasy Points (Season) by Player Age: Fullbacks"
  ) +
  theme_classic() +
  theme(legend.position = "none")

ggplotly(
  plot_rawFantasyPointsByAgeFB,
  tooltip = c("age","fantasy_points","text","label"))

Figure 11.7: Plot of Raw Trajectories of Fantasy Points by Age for Fullbacks.

11.3.2.3 Running Backs

A plot of Running Backs’ raw fantasy points data by age is in Figure 11.8.

plot_rawFantasyPointsByAgeRB <- ggplot(
  data = player_stats_seasonal_offense_subset %>% 
    filter(position == "RB") %>% 
    mutate(
      age = round(age, 2),
      fantasy_points = round(fantasy_points, 2)
    ),
  mapping = aes(
    x = age,
    y = fantasy_points,
    color = player_id)) +
  geom_line(
    aes(
      text = player_display_name, # add player name for mouse over tooltip
      label = season # add season for mouse over tooltip
  )) +
  scale_color_viridis(discrete = TRUE) +
  labs(
    x = "Player Age (years)",
    y = "Fantasy Points (Season)",
    title = "Fantasy Points (Season) by Player Age: Running Backs"
  ) +
  theme_classic() +
  theme(legend.position = "none")

ggplotly(
  plot_rawFantasyPointsByAgeRB,
  tooltip = c("age","fantasy_points","text","label"))

Figure 11.8: Plot of Raw Trajectories of Fantasy Points by Age for Running Backs.

11.3.2.4 Wide Receivers

A plot of Wide Receivers’ raw fantasy points data by age is in Figure 11.9.

plot_rawFantasyPointsByAgeWR <- ggplot(
  data = player_stats_seasonal_offense_subset %>% 
    filter(position == "WR") %>% 
    mutate(
      age = round(age, 2),
      fantasy_points = round(fantasy_points, 2)
    ),
  mapping = aes(
    x = age,
    y = fantasy_points,
    color = player_id)) +
  geom_line(
    aes(
      text = player_display_name, # add player name for mouse over tooltip
      label = season # add season for mouse over tooltip
  )) +
  scale_color_viridis(discrete = TRUE) +
  labs(
    x = "Player Age (years)",
    y = "Fantasy Points (Season)",
    title = "Fantasy Points (Season) by Player Age: Wide Receivers"
  ) +
  theme_classic() +
  theme(legend.position = "none")

ggplotly(
  plot_rawFantasyPointsByAgeWR,
  tooltip = c("age","fantasy_points","text","label"))

Figure 11.9: Plot of Raw Trajectories of Fantasy Points by Age for Wide Receivers.

11.3.2.5 Tight Ends

A plot of Tight Ends’ raw fantasy points data by age is in Figure 11.10.

plot_rawFantasyPointsByAgeTE <- ggplot(
  data = player_stats_seasonal_offense_subset %>% 
    filter(position == "TE") %>% 
    mutate(
      age = round(age, 2),
      fantasy_points = round(fantasy_points, 2)
    ),
  mapping = aes(
    x = age,
    y = fantasy_points,
    color = player_id)) +
  geom_line(
    aes(
      text = player_display_name, # add player name for mouse over tooltip
      label = season # add season for mouse over tooltip
  )) +
  scale_color_viridis(discrete = TRUE) +
  labs(
    x = "Player Age (years)",
    y = "Fantasy Points (Season)",
    title = "Fantasy Points (Season) by Player Age: Tight Ends"
  ) +
  theme_classic() +
  theme(legend.position = "none")

ggplotly(
  plot_rawFantasyPointsByAgeTE,
  tooltip = c("age","fantasy_points","text","label"))

Figure 11.10: Plot of Raw Trajectories of Fantasy Points by Age for Tight Ends.

11.3.3 Linear Regression Models

11.3.3.1 Null Model

pointsPerSeason_nullModel <- lm(
  fantasy_points ~ 1,
  data = player_stats_seasonal_offense_subset,
  na.action = "na.exclude"
)

summary(pointsPerSeason_nullModel)

Call:
lm(formula = fantasy_points ~ 1, data = player_stats_seasonal_offense_subset, 
    na.action = "na.exclude")

Residuals:
   Min     1Q Median     3Q    Max 
-68.13 -53.35 -29.35  28.75 418.61 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  60.8540     0.6321   96.27   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 73.53 on 13531 degrees of freedom
  (1043 observations deleted due to missingness)

summary(pointsPerSeason_nullModel)$r.squared

[1] 0

AIC(pointsPerSeason_nullModel)

[1] 154719.7

MuMIn::AICc(pointsPerSeason_nullModel)

[1] 154719.7

A plot of the model-implied trajectories of fantasy points by age from the null model is in Figure 11.11.

pointsPerSeason_positionAge_newData$fantasyPoints_nullModel <- predict(
  object = pointsPerSeason_nullModel,
  newdata = pointsPerSeason_positionAge_newData
)

ggplot2::ggplot(
  data = pointsPerSeason_positionAge_newData,
  mapping = aes(
    x = age,
    y = fantasyPoints_nullModel
  )
) + 
  geom_line(linewidth = 2) +
  labs(
    x = "Player Age (years)",
    y = "Fantasy Points (Season)",
    title = "Fantasy Points (Season) by Player Age and Position",
    subtitle = "Null Model"
  ) +
  theme_classic()

Figure 11.11: Plot of Model-Implied Trajectories of Fantasy Points by Age in Null Model.

11.3.3.2 Linear Model

pointsPerSeason_linearRegression <- lm(
  fantasy_points ~ positionFactor + ageCentered20 + positionFactor:ageCentered20,
  data = player_stats_seasonal_offense_subset,
  na.action = "na.exclude"
)

summary(pointsPerSeason_linearRegression)

Call:
lm(formula = fantasy_points ~ positionFactor + ageCentered20 + 
    positionFactor:ageCentered20, data = player_stats_seasonal_offense_subset, 
    na.action = "na.exclude")

Residuals:
    Min      1Q  Median      3Q     Max 
-163.81  -51.22  -16.91   39.68  357.82 

Coefficients:
                               Estimate Std. Error t value Pr(>|t|)    
(Intercept)                     11.0884    14.3594   0.772  0.44002    
positionFactorQB                95.0782    15.2244   6.245 4.51e-10 ***
positionFactorRB                68.3543    15.0568   4.540 5.73e-06 ***
positionFactorTE                14.8749    15.1809   0.980  0.32720    
positionFactorWR                42.1981    14.8233   2.847  0.00443 ** 
ageCentered20                    0.7358     1.9891   0.370  0.71144    
positionFactorQB:ageCentered20   2.1806     2.0611   1.058  0.29012    
positionFactorRB:ageCentered20  -1.1383     2.1181  -0.537  0.59102    
positionFactorTE:ageCentered20   1.3905     2.1044   0.661  0.50879    
positionFactorWR:ageCentered20   1.1822     2.0663   0.572  0.56726    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 75.72 on 6435 degrees of freedom
  (8130 observations deleted due to missingness)
Multiple R-squared:  0.1423,    Adjusted R-squared:  0.1411 
F-statistic: 118.6 on 9 and 6435 DF,  p-value: < 2.2e-16

summary(pointsPerSeason_linearRegression)$r.squared

[1] 0.1422979

AIC(pointsPerSeason_linearRegression)

[1] 74077.77

MuMIn::AICc(pointsPerSeason_linearRegression)

[1] 74077.81

A plot of the model-implied trajectories of fantasy points by age from the linear regression model is in Figure 11.12.

pointsPerSeason_positionAge_newData$fantasyPoints_linearRegression <- predict(
  object = pointsPerSeason_linearRegression,
  newdata = pointsPerSeason_positionAge_newData
)

ggplot2::ggplot(
  data = pointsPerSeason_positionAge_newData,
  mapping = aes(
    x = age,
    y = fantasyPoints_linearRegression,
    color = positionFactor
  )
) + 
  geom_line(linewidth = 2) +
  labs(
    x = "Player Age (years)",
    y = "Fantasy Points (Season)",
    title = "Fantasy Points (Season) by Player Age and Position",
    subtitle = "Linear Regression Model",
    color = "Position"
  ) +
  theme_classic()

Figure 11.12: Plot of Model-Implied Trajectories of Fantasy Points by Age in Linear Regression Model.

11.3.3.3 Quadratic Model

pointsPerSeason_quadraticRegression <- lm(
  fantasy_points ~ positionFactor + ageCentered20 + ageCentered20Quadratic + positionFactor:ageCentered20 + positionFactor:ageCentered20Quadratic,
  data = player_stats_seasonal_offense_subset,
  na.action = "na.exclude"
)

summary(pointsPerSeason_quadraticRegression)

Call:
lm(formula = fantasy_points ~ positionFactor + ageCentered20 + 
    ageCentered20Quadratic + positionFactor:ageCentered20 + positionFactor:ageCentered20Quadratic, 
    data = player_stats_seasonal_offense_subset, na.action = "na.exclude")

Residuals:
    Min      1Q  Median      3Q     Max 
-197.10  -50.92  -17.12   39.78  357.46 

Coefficients:
                                        Estimate Std. Error t value Pr(>|t|)
(Intercept)                              -3.2343    33.8844  -0.095   0.9240
positionFactorQB                        145.6404    35.1887   4.139 3.54e-05
positionFactorRB                         83.7797    34.9318   2.398   0.0165
positionFactorTE                         25.8075    35.1789   0.734   0.4632
positionFactorWR                         51.6417    34.6205   1.492   0.1358
ageCentered20                             5.1428     9.6526   0.533   0.5942
ageCentered20Quadratic                   -0.2954     0.6331  -0.467   0.6408
positionFactorQB:ageCentered20          -11.4153     9.8800  -1.155   0.2480
positionFactorRB:ageCentered20           -5.9440    10.0226  -0.593   0.5532
positionFactorTE:ageCentered20           -1.9854     9.9837  -0.199   0.8424
positionFactorWR:ageCentered20           -1.5417     9.8933  -0.156   0.8762
positionFactorQB:ageCentered20Quadratic   0.7527     0.6412   1.174   0.2404
positionFactorRB:ageCentered20Quadratic   0.3247     0.6613   0.491   0.6235
positionFactorTE:ageCentered20Quadratic   0.2308     0.6515   0.354   0.7232
positionFactorWR:ageCentered20Quadratic   0.1767     0.6501   0.272   0.7858
                                           
(Intercept)                                
positionFactorQB                        ***
positionFactorRB                        *  
positionFactorTE                           
positionFactorWR                           
ageCentered20                              
ageCentered20Quadratic                     
positionFactorQB:ageCentered20             
positionFactorRB:ageCentered20             
positionFactorTE:ageCentered20             
positionFactorWR:ageCentered20             
positionFactorQB:ageCentered20Quadratic    
positionFactorRB:ageCentered20Quadratic    
positionFactorTE:ageCentered20Quadratic    
positionFactorWR:ageCentered20Quadratic    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 75.62 on 6430 degrees of freedom
  (8130 observations deleted due to missingness)
Multiple R-squared:  0.1451,    Adjusted R-squared:  0.1433 
F-statistic: 77.98 on 14 and 6430 DF,  p-value: < 2.2e-16

summary(pointsPerSeason_quadraticRegression)$r.squared

[1] 0.1451436

AIC(pointsPerSeason_quadraticRegression)

[1] 74066.35

MuMIn::AICc(pointsPerSeason_quadraticRegression)

[1] 74066.44

A plot of the model-implied trajectories of fantasy points by age from the regression model with a quadratic term for age is in Figure 11.13.

pointsPerSeason_positionAge_newData$fantasyPoints_quadraticRegression <- predict(
  object = pointsPerSeason_quadraticRegression,
  newdata = pointsPerSeason_positionAge_newData
)

ggplot2::ggplot(
  data = pointsPerSeason_positionAge_newData,
  mapping = aes(
    x = age,
    y = fantasyPoints_quadraticRegression,
    color = positionFactor
  )
) + 
  geom_line(linewidth = 2) +
  labs(
    x = "Player Age (years)",
    y = "Fantasy Points (Season)",
    title = "Fantasy Points (Season) by Player Age and Position",
    subtitle = "Quadratic Regression Model",
    color = "Position"
  ) +
  theme_classic()

Figure 11.13: Plot of Model-Implied Trajectories of Fantasy Points by Age in Quadratic Regression Model.

11.3.3.4 Compare Models

anova(
  pointsPerSeason_linearRegression,
  pointsPerSeason_quadraticRegression
)

AIC(
  pointsPerSeason_nullModel,
  pointsPerSeason_linearRegression,
  pointsPerSeason_quadraticRegression
  )

lmModels <- list(
  "nullModel" = pointsPerSeason_nullModel,
  "linearRegression" = pointsPerSeason_linearRegression,
  "quadraticRegression" = pointsPerSeason_quadraticRegression
)

bbmle::AICtab(lmModels)

                    dAIC    df
quadraticRegression     0.0 16
linearRegression       11.4 11
nullModel           80653.3 2 

summary(pointsPerSeason_nullModel)$r.squared

[1] 0

summary(pointsPerSeason_linearRegression)$r.squared

[1] 0.1422979

summary(pointsPerSeason_quadraticRegression)$r.squared

[1] 0.1451436

deviance(pointsPerSeason_nullModel)

[1] 73167060

deviance(pointsPerSeason_linearRegression)

[1] 36895690

deviance(pointsPerSeason_quadraticRegression)

[1] 36773276

logLik(pointsPerSeason_nullModel)

'log Lik.' -77357.85 (df=2)

logLik(pointsPerSeason_linearRegression)

'log Lik.' -37027.89 (df=11)

logLik(pointsPerSeason_quadraticRegression)

'log Lik.' -37017.18 (df=16)

11.3.4 Mixed Models

By accounting for which player each observation comes from using mixed models, we can examine the association between age and fantasy points in a more meaningful way, without violating the assumption in multiple regression that the observations are independent (i.e., that the residuals are uncorrelated).

11.3.4.1 Random Intercepts Model

pointsPerSeason_randomIntercepts <- lmerTest::lmer(
  fantasy_points ~ 1 + (1 | player_idFactor),
  data = player_stats_seasonal_offense_subset,
  REML = FALSE,
  control = lmerControl(optimizer = "bobyqa")
)

summary(pointsPerSeason_randomIntercepts)

Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
  method [lmerModLmerTest]
Formula: fantasy_points ~ 1 + (1 | player_idFactor)
   Data: player_stats_seasonal_offense_subset
Control: lmerControl(optimizer = "bobyqa")

     AIC      BIC   logLik deviance df.resid 
148044.0 148066.6 -74019.0 148038.0    13529 

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-6.5321 -0.4328 -0.1682  0.3538  5.5810 

Random effects:
 Groups          Name        Variance Std.Dev.
 player_idFactor (Intercept) 2197     46.87   
 Residual                    2335     48.32   
Number of obs: 13532, groups:  player_idFactor, 3358

Fixed effects:
             Estimate Std. Error        df t value Pr(>|t|)    
(Intercept)   44.6942     0.9595 3834.9436   46.58   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

MuMIn::r.squaredGLMM(pointsPerSeason_randomIntercepts)

     R2m       R2c
[1,]   0 0.4847453

performance::icc(pointsPerSeason_randomIntercepts)

AIC(pointsPerSeason_randomIntercepts)

[1] 148044

AICc(pointsPerSeason_randomIntercepts)

numeric(0)

A plot of the model-implied trajectories of fantasy points by age from the mixed model with random intercepts is in Figure 11.14.

pointsPerSeason_positionAge_newData$fantasyPoints_randomIntercepts <- predict(
  object = pointsPerSeason_randomIntercepts,
  newdata = pointsPerSeason_positionAge_newData,
  re.form = NA
)

ggplot2::ggplot(
  data = pointsPerSeason_positionAge_newData,
  mapping = aes(
    x = age,
    y = fantasyPoints_randomIntercepts
  )
) + 
  geom_line(linewidth = 2) +
  labs(
    x = "Player Age (years)",
    y = "Fantasy Points (Season)",
    title = "Fantasy Points (Season) by Player Age",
    subtitle = "Random Intercepts Model"
  ) +
  theme_classic()

Figure 11.14: Plot of Model-Implied Trajectories of Fantasy Points by Age in Random Intercepts Mixed Model.

A plot of individuals’ model-implied trajectories of fantasy points by age from the mixed model with random intercepts is in Figure 11.15.

player_stats_seasonal_offense_subsetCC$fantasyPoints_randomIntercepts <- predict(
  object = pointsPerSeason_randomIntercepts,
  newdata = player_stats_seasonal_offense_subsetCC
)

plot_individualFantasyPointsRandomIntercepts <- ggplot(
  data = player_stats_seasonal_offense_subsetCC,
  mapping = aes(
    x = age,
    y = fantasyPoints_randomIntercepts,
    #color = positionFactor,
    group = player_id)) +
  geom_smooth(
    aes(
      x = age,
      y = fantasyPoints_randomIntercepts,
      text = player_display_name, # add player name for mouse over tooltip
      label = season # add season for mouse over tooltip
    ),
    se = FALSE,
    linewidth = 0.5,
    color = "black") +
  geom_line(
    mapping = aes(
      x = age,
      y = fantasyPoints_randomIntercepts
    ),
    data = pointsPerSeason_positionAge_newData,
    inherit.aes = FALSE,
    se = TRUE,
    color = "#3366FF",
    linewidth = 2
  ) +
  geom_point(
    aes(
      x = age,
      y = fantasyPoints_randomIntercepts,
      text = player_display_name, # add player name for mouse over tooltip
      label = season # add season for mouse over tooltip
    ),
    size = 1,
    color = "transparent" # make points invisible but keep tooltips
  ) +
  labs(
    x = "Player Age (years)",
    y = "Fantasy Points (Season)",
    title = "Fantasy Points (Season) by Player Age and Position: Random Intercepts Model",
    #color = "Position"
  ) +
  theme_classic()

ggplotly(
  plot_individualFantasyPointsRandomIntercepts,
  tooltip = c("age","fantasyPoints_randomIntercepts","text","label")
)